function [diff, eq] = solve_cL(D,param,glob,options)

param.D = D;

%% Initialise guess
c0   = zeros(2*glob.Ns,1);

c1old   = reshape(c0(1:glob.Ns), glob.n(1), glob.n(2));
c2old   = reshape(c0(glob.Ns+1:end), glob.n(1), glob.n(2));

cold    = [c1old; c2old];

if isempty(options.cresult)==0
    cold = options.cresult;
end

%% Bellman iterations
%fprintf('~~~~~ Bellman iterations ~~~~~\n');
dc = 10;
counter = 0;
while dc > options.tolc && counter < options.Nbell
    counter = counter + 1;
         
    v     = solve_valfunc(cold,glob.s,param,glob,options);
 
    % begin HPI
   v2     = reshape(v.v2, glob.n(1), glob.n(2));
   v2     = griddedInterpolant(glob.B, glob.S,v2, 'linear');   
   for i = 1:options.hpisteps
       v1     = valfunc_I(v2,glob.s,v.y,param,glob,options);
       v2     = glob.PnK*v1;
       v2     = reshape(v2, glob.n(1), glob.n(2));
       v2     = griddedInterpolant(glob.B, glob.S,v2, 'linear');   
   end
    
   v2 = reshape(glob.PnK*v1, glob.n(1), glob.n(2));
    % end

    v2 = reshape(v2, glob.n(1), glob.n(2));
    v1 = reshape(v1, glob.n(1), glob.n(2));
    c = [v1;v2];
          
%    dc = max(max(max(abs((c - cold)./cold))))
    dc = max(max(max(abs(c-cold))));
    cold = options.damp*cold + (1-options.damp)*c;

end

if dc < options.tolc
   eq.convflag = 1; 
%   disp('~~ bellman iteration converged ~~')
else
   eq.convflag = 0;
   disp('<strong>BELLMAN ITERATION DID NOT CONVERGE</strong>')
   dc
end

%% compute value of entry
v1     = griddedInterpolant(glob.B, glob.S,v1, 'linear');   

ve.v1   = v1(zeros(size(glob.agrid)), glob.agrid); % initial value
eq.ve_l = 0*ve.v1;

eq.ve = ve;

eq.eve  = sum(eq.ve.v1.*glob.a_stat_E');


%% Solve optimal decision again on a finer grid
tmp      = glob.PnK;
glob.PnK = glob.PnKf;
vf       = solve_valfunc(cold,glob.sf,param,glob,options);
glob.PnK = tmp;

%% Get labor demand decision
eq.l = menufun('L',glob.sf,vf.y,param,glob,options);

%% Pack up output
%eq.bp     = bp;
eq.p      = v.p;
eq.y      = v.y;
eq.pf     = vf.p;
eq.yf     = vf.y;
eq.p_exit = vf.p_exit;
%eq.y      = y;
eq.c      = c;
eq.v      = v;
eq.v2     = v2;
%eq.L      = L;
%eq.Q      = Q;
%eq.Lp     = Lp;
%eq.La     = La;
eq.vf     = vf;

diff = eq.eve - param.ce;

end

